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ABSTRACT 

Lagrangian smoothed particle hydrodynamics (SPH) is a well-established approach to 
model fluids in astrophysical problems, thanks to its geometric flexibility and ability 
to automatically adjust the spatial resolution to the clumping of matter. However, a 
number of recent studies have emphasized inaccuracies of SPH in the treatment of 
fluid instabilities. The origin of these numerical problems can be traced back to spuri- 
ous surface effects across contact discontinuities, and to SPH's inherent prevention of 
mixing at the particle level. We here investigate a new fluid particle model where the 
density estimate is carried out with the help of an auxiliary mesh constructed as the 
Voronoi tessellation of the simulation particles instead of an adaptive smoothing ker- 
nel. This Voronoi-based approach improves the ability of the scheme to represent sharp 
contact discontinuities. We show that this eliminates spurious surface tension effects 
present in SPH and that play a role in suppressing certain fluid instabilities. We find 
that the new 'Voronoi Particle Hydrodynamics' described here produces comparable 
results than SPH in shocks, and better ones in turbulent regimes of pure hydro- 
dynamical simulations. We also discuss formulations of the artificial viscosity needed 
in this scheme and how judiciously chosen correction forces can be derived in order 
to maintain a high degree of particle order and hence a regular Voronoi mesh. This 
is especially helpful in simulating self-gravitating fluids with existing gravity solvers 
used for N-body simulations. 
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1 INTRODUCTION 

Numerical simulations have become an important research 
tool in many areas of astrophysics, in particular in cosmic 
structure formation and galaxy formation. This is in part 
because the physical conditions involved cannot be repro- 
duced in laboratories on Earth, so that simulations serve as 
a replacement for experiments. Perhaps more importantly, 
simulations in principle allow a full modeling of all the in- 
volved physics. However, a significant problem in practice is 
that that the equations one wants to solve first have to be 
numerically discretized in a suitable fashion. The accuracy 
of simulations depends strongly on the properties of this dis- 
cretization, and it hence remains an important task to find 
improved numerical schemes for astrophysical applications. 

In cosmic structure formation, matter is initially essen- 
tially uniformly distributed, but clusters with time under 
the action of self-gravity to enormous density contrasts, pro- 
ducing galaxies of vastly different sizes. Given the variety of 
involved geometries, densities and velocities, it is clear that 
a Lagrangian method, where the mass of a resolution ele- 
ment stays (roughly) constant, would be most convenient. 
This is because a Lagrangian method automatically concen- 
trates the resolution in regions where the galaxies form, and 



hence focuses the numerical effort on the regions of interest. 
On the other hand, traditional mesh-based approaches to 
hydrodynamics, so-called Eulerian methods, discretize the 
volume in a set of cells and do not follow the clustering of 
matter, unless this is attempted with a suitable adaptive 
mesh-refinement strategy. 

The by far most widely used Lagrangian approach 
in structure formation is smoothed partic l e hydrodynam- 
ics (S PH, as reviewed by iMonagharJl 19921 . 2005; Ro sswod 
2009), a technique that dates back to particle-based ap- 
proa ches first developed in astronomy more than 30 year s 
ago (|Lucvlll977l ; iGingold fc Monaghanlll977l ; lLarsonlll978h . 
In this method the fluid is discretized in terms of par- 
ticles of fixed mass, which are used to construct an ap- 
proximation to the Euler equations based on the adap- 
tive kernel interpolation technique. SPH can be very easily 
coupled to self-gravity, it is remarkably robust (e.g. neg- 
ative densities cannot arise), and the introduction of extra 
physics (e.g. feedback processes in the context of star forma- 
tion) is intuitive. All of these properties have made it very 
popular for p roblems such as planet formation or galaxy 
mergers (e.g. iMihos fc Hernquistl Il99rj; iMaver et al.l [2003 ; 



iRobertson et all |2004| ; iDolag et al.l l2005n ~ where spatially 
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separated regions of the simulation volume feature widely 
different densities. 

However, recent studies have highlighted a number of 
differences in the results of SPH-based calculations com- 
pared to more traditional grid-based Eulerian methods for 
hydrodynamics. For example, the two methods appear to 
disagree about the entropy produced in the central region of 
a forming galaxy cluster under non-radiative conditions, as 
first seen in the 'S anta Barbara cluster comparison project' 
l|Frenk et al.lll999h . It has been suggested that this problem 
may be caused by a suppression of the Raleigh- Taylor fluid 
instability in SPH (IMitchell et al.ll2009r ) and the lack of mix- 



ingat the part icle level dTasker et aljboosl ; IWadslev et al.l 
|200ST ). Indeed, lAgertz et all 1120071 ) have shown that SPH 
tends to suppress Kelvin-Helmholtz fluid instabilities in 
shear flows across interfaces with sizable density jumps. In 
such a situation, SPH's density estimate leads to spurious 
forces at the interface which produce an artificial 'gap' in the 
particle distribution and a surface tension effect that ulti- 
mately produces errors in the hydrodynamical evolution. To 
what extent these numerical artifacts negatively affect the 
global accuracy of simulations in practice is unclear, and this 
can in any case be expected to be problem dependent. How- 
ever, an improvement of standard SPH that avoids these 
errors is obviously desirable. 

F irst p r oposa ls in this direction have recently been 
made. IPricel (|2008l ) suggests to introduce artificial heat con- 
duction into SPH such that discontinuities in the tempera- 
ture field are smoothed out, in analogy to the ordinary artifi- 
cial viscosity that effectively smoothes out discontinuities in 
the velocity field occurring at shocks. This heat conduction 
produces a soft instead of an abrupt transition of the specific 
entropy across a contact discontinuity, which in turn helps 
to better represent the growth of Kelvi n-Helmholtz insta- 
bilities at such interfaces. More recently. | Read et aL (|2009T l 
have modified an idea by iRitchie fc Thomas! ifeooih for a 
modified SPH density estimate that assumes that the local 
neighbours have similar pressures, and which is designed to 
avoid the 'pressure blip' in the standard approach at con- 
tact discontinuities. Together with a modified kernel shape 
and a drastically enlarged number of neighbours (by a fac- 
tor of ~ 10, imply i ng a s imilar increase in the computational 
cost). iRead et all (|20 09 : 1 obtained better growth of Kelvin- 
Helmholtz instabilities across density jumps. 

In this work we follow a different approach that elimi- 
nates the ordinary SPH kernel altogether. Instead, we use 
the distribution of points with variable masses to construct 
an auxiliary mesh, which is then used to derive local density 
estimates. If the particle hydrodynamics is derived from a 
Lagrangian, it turns out that obtaining this density estimate 
is already sufficient to uniquely determine the equations 
of motion. The use of Delaunay tessellations to construct 
density fields from ar bitrary point sets has been dis- 
cussed in the literature l|Schaap fc van de Wevgaertll2000al; 



Icke fc van de Wevgaert Il987l ; van de Weygaert fc Ickd 
19891 ; iPelupessv et al l l2003h , but as we show in this 



paper, its topological dual, the Voronoi tessellation, is 
actually preferable for our hydrodynamical application. In 
the Voronoi tessellation, to every particle a polyhedra is 
assigned which encompasses the space closer to this particle 
than to any other. Based on these volumes associated with 
each particle, local densities and hydrodynamical forces 



can be estimated, leading to an interesting alternative to 
SPH. In particular, it is immediately clear that unlike SPH 
this approach yields a consistent discretization not only 
of the mass but also of the volume, which should help to 
yield an improved representation of contact discontinu- 
ities. We note that a conceptionally similar approach to 
Vo ronoi based particl e hydr odynamics was first discussed 
bv lSerrano fc Espanoll (|200ll ) in the context of a mesoscopic 
fluid particle model. We here extend this idea to the 
treatment of the Euler equations in astrophysical systems. 

We emphasize that the method we introduce in this 
study is radically diffe rent from to t he one implemented in 
the new AREPO code l|Springelll2009l 'l. Whereas the latter is 
also based on a (moving) Voronoi tessellation, it employs a 
finite volume scheme with a Riemann solver to compute hy- 
drodynamical fluxes across mesh boundaries. This involves 
an explicit second-order reconstruction of the fluid through- 
out the volume, and allows for changes of the mass contained 
in each cell even if the mesh is on average moving with the 
flow. In contrast, we here derive a fluid particle model from 
a discretized Lagrangian in which the masses of each ele- 
ment stay strictly constant, and in which the motion of the 
particles is governed by pairwise pressure force exchanged 
between them. While AREPO is conceptually close to the 
techniques used in Eulerian hydrodynamics, the method we 
study here is conceptually close to SPH. 

This paper is structured as follows. In Section [2] we 
discuss how the equations of motion can be derived for the 
Lagrangian particle approach to hydrodynamics discussed 
here. We will also present suitable formulations of artificial 
viscosity for our scheme. In Section we discuss the role of 
the regularity of Voronoi cells and means to improve it. We 
briefly describe the implementation of our numerical scheme 
in a modified version of the GADGET code in Section 3J and 
then turn in Section[5]to a description of results for a suite of 
test problems with our new 'Voronoi Particle Hydrodynam- 
ics' (VPH) scheme. These tests range from simple shock- 
tube problems, to fluid instabilities, and three-dimensional 
stripping of gas in a supersonic flow. Finally, we summarize 
our conclusions in Section [fj] In two Appendices, we discuss 
gradient operators for Voronoi meshes and give the deriva- 
tion of correction forces that can be used to maintain very 
regular mesh geometries, if desired. 



2 PARTICLE BASED HYDRODYNAMICS 

We begin by introducing our methodology for a particle- 
based fluid dynamics based on Voronoi tessellations. This 
method is close in spirit to SPH, but differs in important 
aspects. Where appropriate, we discuss these differences in 
detail. 



2.1 A Lagrangian approach for particle based 
fluid dynamics 

We discretize the fluid in terms of N mass elements of mass 
rrii. The discretized fluid Lagrangian can then be adopted 
as 



[i 



rm Ui(pi, Si) 



(1) 
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This is simply the difference of the kinetic and thermal en- 
ergy of the particles. The thermal energy Ui per unit mass 
depends both on the density pi and the specific entropy s» 
of the particle. In this work, we aim to approximate invis- 
cid ideal gases, hence the equation of state (EOS) is that 
of a polytropic gas, where the pressure is Pi = Sip 1 , and 
the entropic function Si (or simply 'entropy' for short, since 
it depends only on the thermodynamic entropy) labels the 
adiabat on which this gas element resides. 

When the specific entropy Si and the mass rrn remain 
constant for a fluid element i, the internal energy m changes 
according to |j = Using this result, we can readily write 
down the Lagrangian equations of motion of (reversible) 
fluid dynamics: 



rrnri 



' dri 

duj dpj 
' dpj dr l 

EPj dpj dVj 
m3 ~p\W^ 



3 
N 



3 
N 



(2) 



We see that the primary input required for a more explicit 
form of the equations is a density estimate based on the par- 
ticle coordinates, or alternatively, an estimate of the volume 
associated with a given particle. 

SPH addresses this task with a kernel estimation tech- 
nique to obtain the density, where an adaptive spherically 
symmetric smoothing kernel is employed to calculate the 
density based on the spatial distribution of an approximately 
fixed number of nearest neighbours. The Lagrangian then 
uniquely determines the equations of motion that simulta- 
neou sly conserve energy and entropy (|Springel fe Hernquistl 
2002). However, we note that SPH does not achieve a con- 
sistent volume estimate, i.e. the sum of the effective volumes 
of the particles, Vi = rrii/pi, is not guaranteed to be equal 
to the total simulated volume. Furthermore, the inherent 
smoothing operation in the density estimate is bound to be 
inaccurate at contact discontinuities and phase interfaces, 
where the density may discontinuously jump by a large fac- 
tor. In the following, we therefore look for alternative ways to 
construct density estimates which improve on these deficits. 



2.2 Density estimates with tessellation techniques 

One promising approach for more accurate density and spe- 
cific volume estimates lies in the use of an auxiliary mesh 
that is generated by the particle distribution. A mesh can 
readily yield a partitioning of the volume such that the to- 
tal volume is conserved, and also allows multiple ways to 
'spread out' the particle masses rrii in a conservative fashion 
such that an estimate of the density field is obtained. 

There are two basic geometric constructions that sug- 
gest thems elves as s uch m esh candidates. These are the 
Delaunay (jPirichletl Il850l ) and the Voronoi tessellations 
|Voronoilll908l )7 which are in fact mathematically closely re- 
lated, as we discuss below. In the Voronoi tessellation, space 
is subdivided into non-overlapping polyhedra which each 
encompass the volume which is closer to its corresponding 



point than to any other point. The surfaces of these poly- 
hedra are therefore the bisectors to the nearest neighbours. 
The Delaunay tessellation on the other hand decomposes 
space into a set of tetrahedra (or triangles in 2D), with ver- 
tices at the point coordinates. The defining property of the 
Delaunay tessellation is that the circumcircles of the tetra- 
hedra do not contain any of the points in their interior. This 
property in fact makes this tessellation uniquely determined 
for points in general position. 

It turns out that these two tessellations are dual to 
each other; to each edge of the Delaunay tessellation cor- 
responds a face of the Voronoi tessellation, and the circum- 
circle centres of the Delaunay tetrahedra are the vertices of 
the Voronoi faces. One implication of this is that the Delau- 
nay and Voronoi Tessellations can be easily transformed into 
each other. In practice it is typically simpler to always con- 
struct the Delaunay tessellation, even if one works with the 
Voronoi, because the former has more efficient and simpler 
algorithms for construction. 

Both tessellati ons can in principle be us e d to d erive 
density estimates. ISchaap fc van de Weygaertl l|2000bh in- 
troduced the D elaunay Tessellat i on Fie ld Estimator (DTFE) 
technique, and lPehxpessv et"ai1 (|2003l ) showed that it offers 
superior resolution compared to SPH-like density estimates 
for detecting cosmological large-scale structure (for a review 
see Ivan de Weygaert fc Sc haap 200jl). In this approach, the 
total volume of the contiguous set of Delaunay cells around 
a point is used to assign particle densities, and a full density 
field can be constructed by linearly interpolating the densi- 
ties inside each Delaunay tetra hedron. As a possible a ppli- 
cation of this density estimate, iPelupessv et al.l (|2003t ) also 
suggested its use in a particle-based hydrodynamic scheme. 
However, we caution that a rather serious short-coming of 
the Delaunay tessellation in this context is that the tessella- 
tion may occasionally change discontinuously as a function 
of the particle coordinates. This happens whenever a parti- 
cle moves over the circumcircle of one of the tetrahedra. An 
infinitesimal particle motion can hence be sufficient to cre- 
ate finite changes in the volume of its associated contiguous 
Delaunay cell (this is the union of all Delaunay tetrahedra 
of which the given point is one of the vertices) of a parti- 
cle. As a result, the thermal energy of the point set is not a 
continuous function of the particle coordinates. This makes 
the DTFE technique ill suited to be the basis of a hydrody- 
namical particle method. 

On the other hand, the volumes of the Voronoi cells 
always depend continuously on the particle coordinates, de- 
spite the fact that topological changes of the tessellation 
may occur as a result of particle motion. This is because 
flips of edges in the Delaunay tessellation happen precisely 
when the corresponding Voronoi faces have vanishing area. 
Another advantage of the Voronoi tessellation is that it is 
always uniquely defined for any distribution of the points, 
whereas for certain degenerate point sets (those where more 
than four points lie on a common circumsphere) , more than 
one valid Delaunay tessellation may exist, which can then 
make Delaunay-derived density estimates non-unique. We 
remark that the uniqueness of the Voronoi tessellation does 
not hold in reverse, i.e. a given Voronoi tessellation can in 
general be produced by a number of different point distri- 
butions. This has important consequences for the stability 
of the scheme, as we discuss later on in more detail. 
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Figure 1. Section of a Voronoi diagram for a set of points marked 
with asterisks. All the triangles in the dual Delaunay tessellation 
that are shared by the point in the centre are marked in blue. 
The vectors cy and ey needed to calculate the derivative of the 
volume are also marked. 



Based on the above, the Voronoi tessellation is a promis- 
ing construction for a particle fluid model, hence we adopt it 
in the following. In particular, we shall associate the volume 
Vi of a Voronoi cell with its corresponding point, yielding a 
consistent decomposition of the total simulated volume. The 
simplest possible density estimate is then simply given by 



(3) 



which we shall use in this paper. More involved higher-order 
density field reconstructions could be considered as well, an 
idea we leave for future work. 

Based on this density estimate, and given the specific 
entropies of each particle, the local pressure and the ther- 
mal pressure per unit mass can be computed. Also, one may 
define a gradient operator for the Voronoi mesh (see our dis- 
cussion in Appendix I Al|) . which could be used to estimate 
pressure gradients, and hence to yield discretizations of the 
Euler equations. However, a better approach is to start from 
the discretized Lagrangian, as this automatically gives equa- 
tions of motion that satisfy the conservation laws. We shall 
adopt this strategy in the following. 



2.3 Equations of motion for Voronoi-based 
particle hydrodynamics 

Since the volumes of the Voronoi cells depend only on the 
configuration of the points, we can readily obtain the equa- 
tions of motion if we find the partial derivative of a cell 
volume with respect to an y of the particle coo r dinat es. We 
here adopt the result of ISerrano fc Esoanoll (|200lh . who 
showed that the releva nt derivative is given by (see also 
|Pe Fabritiis et alj|2006h 



dri 



-Ai 



for i / j, 



(4) 



where denotes the gradient operator with respect to Vi. 
Here Rij is the distance between two neighboring points, 
e ij = ( r i ~ r i)/Rij denotes a unit vector from i to the 
neighbour j, which is normal to the Voronoi face of area Aij 
between cells i and j. Formally, we can define Aij for any pair 
of different particles, but if i and j are not neighbours in the 
Voronoi tessellation (i.e. do not share a face), we set Aij = 0, 
implying that in sums that involve the factor Ay only the 
direct neighbours contribute. Note that equation @ holds 
only for j i. But one can readily derive an expression for 
dVi/dri by invoking volume conservation. This yields 



dVj 
dri 



(5) 



As sketched in Figure [TJ the vector Cy points from the mid- 
point between i and j to the centroid of the face Aij, and 
is orthogonal to eij. The term involving ey/2 can be easily 
understood geometrically from the change of the volumes of 
the pair of pyramids spanned by the face between i and j 
and the two points. But if the center of the face is displaced 
from the line connecting i and j, a second term involving cy 
appears that stems from the turning of the face when the 
points are moved. 

We are now in a position to write down the resulting 
equations of motion, based on equations @, (JU) and ©. 
This first yields 



dpj 
dri 



VOj_ 



(6) 



which then gives rise to the equations of motion in the form 



(7) 



m t ri = \_. Aij(P, - Pj) 

This is a rather intuitive result, as it shows that motions are 
generated by the pressure differences that occur across faces 
of the tessellation. If the pressures are all equal, the forces 
vanish exactly, unlike in ordinary SPH. 

In the form of equation @, it is not obvious whether the 
forces between a given pair of particles are antisymmetric. 
However, noting the identity ^2j=u Aij£ij = 0, which follows 
from Gauss' theorem, we can restore manifest antisymme- 
try in the equations of motion, which is in general prefer- 
able for numerical reasons. To this end, we simply subtract 
Pi Y\ ,, Aij eij = from (I7J), yielding our final equations of 
motion as 



rrnri 



Ri 



(8) 



which is now pairwise antisymmetric. Note that whereas for- 
mally the sums appearing in these equations are carried out 
over all particles, only the direct neighbours actually con- 
tribute, and these are known from the tessellation. In fact, 
the list of interacting particle pairs is exactly given by the list 
of edges of the underlying Delaunay tessellation, or equiva- 
lently, by the list of faces of the Voronoi tessellation. 

We further note that since the equations of motion have 
been derived from the Lagrangian given in equation (fl]), 
these equations conserve energy, momentum and entropy 
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exactly. In the present form they are hence a description 
of the reversible, adiabatic parts of a flow, but they do not 
yet contain any dissipation, which is however needed to treat 
shocks. If no such dissipation is included, shocks will lead to 
unphysical ringing and oscillations in the fluid. 

2.4 Artificial viscosity 

We follow the s tanda r d SPH ap proach (e.g. 
iGingold fc Monaghanl 1 19771 ; iBalsaral 1 19951 ) and invoke 
an artificial dissipation in the form of an extra friction force 
that reduces the kinetic energy and transforms it into heat. 
There is great freedom in the form of this viscous force, 
but ideally it should only become active where it is really 
needed, i.e. in shocks, and should be negligible away from 
shocks, such that inviscid behavior is ensured there. The 
most widely used and tested formulation of the viscous 
acceleration in SPH schemes is given by 



pair of interacting particles, i.e. (P v i s . 



(dvisc)i 



= ^ vu\ 



Pij 



fi + h 



fi 



IV ■ v\ 



\V-v\i + IV x vU +e 



(9) 
(10) 
(11) 
, (12) 



provided that Vit ■ < 0, i.e. the neighboring particles ap- 
proach each other, otherwise the viscous force that is medi- 
ated by the viscous tensor 11^ is set to zero. In this notation, 
qij represents the difference and cjij the average between the 
quantities q associated with particles i and j. The parame- 
ter e is a tiny value introduced to guard against numerical 
divergences. The parameters a and j3 set the strength of the 
viscosity and are typically set to of order ~ 1. The factors fi 
measure the strength of the local velocity dispersion relative 
to the local shear, and are introduced as so-called Balsara 
switch to reduce the vi scosity if the local flow is dominated 
by shear i|Balsaralll995l ). 

The above formulation of a viscous force can be adopted 
to the Voronoi scheme in a number of ways. We first define 
the projected pairwise velocity as 

wij = «^ry > (13) 



and make the replacement [j,ij — > Wij. This is effectively 
yielding the 'sign al velocity' form of the standard viscosity 
(|Monaghanl I1997T ). For simplicity, we shall also adopt the 
common choice /3 = 2a. We next recognize that in SPH the 
viscous tensor is introduced into the equations of motion 
as if it was an extra pressure of the form P v j sc = |pyHij 
jSpringejlioblh . Using this analogy, and guiding ourselves 
by the form of the Voronoi-based equations of motion @, 
we can readily write down a parameterization of the viscous 
force acting on a particle as 



m<i(Q'visc*)i — ^ ^ A 



—2 -pr e ij 



(14) 



Here we have only introduced a viscous force component 
parallel to the line connecting the two particles, since we 
assume that the 'viscous extra pressure' is the same for a 



(-Pvisc)i- A more 
explicit form of the viscous acceleration is given by the fol- 
lowing expression: 



(ivisc )i 



= a y — 

TO; 

3 



ij 



[Wij Cij 



- 2wij 



(15) 



Note that the viscous force is pairwise antisymmetric, and 
will only become active if two particles approach each other. 
We also want to stress that artificial viscosity parameteriza- 
tions different from that of equation (|15[) are of course pos- 
sible. We here simply adopt this form as a first best guess, 
based on the analogy with the widely tested SPH formula- 
tion. 

It is interesting to compare the artificial viscosity with 
the viscosity terms of the Navier Stokes equation, 

Dv 

m— = -VP + 77 Aw + A V(Vv). (16) 

Neglecting the shear viscosity r\ and approximating the gra- 
dient operator with its Voronoi discretized form (see Ap- 
pendix A), this becomes 



i 

Additionally approximating 



(17) 



(18) 



yields a term like the one linear in w in equation (|15[1 , adding 
some further justification to this form of the viscous force, 
which has the form of an artificial bulk viscosity. 

In order to maintain energy conservation, heat must be 
produced at a rate dE/dt that exactly balances the loss of 
kinetic energy due to the extra friction from the artificial vis- 
cosity. We inject this energy symmetrically into the specific 
entropies of the two particles. Defining the pairwise viscous 
forces as 



(/visc)»j — ~AijP 



the heating is given by 
duj 

dt 2 



<y ^ ^ (jf visc)^ 



(19) 



(20) 



With m — Sip] V(7 — 1) this yields for the rate of entropy 
production 

dsi 
It 



7- 



1 \ " 



if, 



(21) 



In this form, the equations still conserve total energy and 
momentum, while the change of the total entropy is positive 
definite. 

The artificial viscosity is necessary to capture shocks 
and to damp postshock oscillations in the vicinity of shocks, 
but everywhere else in the fluid it can induce spurious dissi- 
pation that distorts the physics of an inviscid gas. In order 
to reduce the influence of the viscosity in regions away from 
shocks, the prefactor a tha t sets the strength of the viscosit y 
can be chosen adaptively (|Morrisl [l997l ; iDolag et all 120051 ) . 
The idea of this dynamic viscosity is that every particle 
gets an individual viscosity strength a which is evolved in 
time according to the differential equation 
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da 
~dt 



+ s. 



(22) 



Here a is decaying to a minimum a* on a timescale r, and 
is increased by the source term S. One possible choice for 
this source term is 



s = h C v • v, 



(23) 



which we adopt in our implementation of a time-variable 
artificial viscosity, using the discretized estimate of the di- 
vergence described in Appendix A. Here both the response 
coefficient £ and the timescale r have to be calibrated em- 
pirically. When a shock arrives in an unperturbed area, a is 
at its minimum and needs to jump very quickly to a higher 
level in order to capture the shock and prevent post shock 
oscillations, whereas behind the shock, the viscosity should 
quickly return to a low value. However, a too large value for 
£ may trigger high viscosity due to the often noisy estimates 
of the V • v term, and if r is too small, the viscosity may 
decay too quickly to capture the shock properly. Finally, the 
minimum viscosity a* can be set to a non-zero value to im- 
prove particle order and thereby reduce noise, at the cost of 
introducing some minimum viscosity. 



Here k = 2n/L is the wavenumber of the Kelvin-Helmholtz 
mode. We shall assume that L is of order the cell dimen- 
sion, which is in turn of order the particle separation, i.e. we 
will set L = Rij when a particle pair of separation Rij is 
considered. Similarly, the relevant velocity jump is simply 
the velocity difference projected onto the face between two 
particles, which is normal to their separation vector, hence 



(25) 



We further assume that the fluid mixing on scales below the 
particle cells can be approximately described as a diffusion 
process, operating with diffusion constant D = x^/^kh, 
where x is a dimensional efficiency that controls the strength 
of the mixing (and which needs to be determined empir- 
ically). We hence effectively model the mixing with heat 
diffusion of the form du/dt = _DV 2 it. 

U sing the SPH discreti zation of thermal conduction as a 
guide (Jubclgas ct al. 2003), we can readily find a discretiza- 
tion of the heat diffusion for the Voronoi particle discretiza- 
tion. We obtain 
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2.5 Treatment of mixing 

Hydrodynamic simulations are able to follow the advection 
of fluids only down to the resolution scale. But especially 
Lagrangian schemes do not include mixing processes of the 
fluid on sub-resolution scales. In Eulerian codes, such mixing 
is implicit whenever a new averaged thermodynamic state 
for a cell is computed after fluxes of gas have entered or 
left it. This mixing keeps the total energy fixed, but will in 
general raise the entropy of the system. In the Lagrangian 
particle approach of SPH and in the Voronoi approach de- 
veloped here, such mixing effects are, however, entirely sup- 
pressed. The specific entropies of neighbouring particles stay 
constant, except when a shock is present. While this reliably 
eliminates unwanted entropy production from advection er- 
rors, it also prevents the proper subresolution production of 
entropy when small-scale fluid instabilities should mix the 
fluid on the resolution scale and produce homogeneous ther- 
modynamic properties. 

We have therefore tried to model this subgrid mix- 
ing with a heuristic model which conjectures that small- 
scale fluid instabilities, if present, equalize the local tem- 
perature field by mixing. This will then smooth out sharp 
contact discontinuities and also tend to equilibrate the 
specific entropies of the cells. S imilar ideas have rec ently 
been discu s sed b y iPricel (120081 1. 1 Wadslev et all (|2008T ) and 
IShen et all (|2009l ) in the context of SPH, but our approach 
differs in detail. In particular, we restrict the averaging to 
shearing layers, and motivate the timescale for mixing di- 
rectly with the growth timescale of the Kelvin-Helmholtz 
instability on the resolution scale. 

The linear theory growth timescale of a perturbation 
across a contact discontinuity with densities pi and pj that 
exhibits a jump in the tangential velocity of size vfj (a shear 
layer) is given by 



Pi + Pj 



2 k v> 



(24) 



Here we introduced a further factor (1 — /y) in a similar 
manner as in (|12[) . This Balsara-like factor is used to re- 
strict the diffusion only to areas where the compression (as 
measured by |V • v\) is clearly negligible compared to the 
shear (as measured by |V x v\). 

Note that this equation preserves the total thermal en- 
ergy, and heat energy only flows from hotter to colder par- 
ticles. The corresponding rates of entropy change for each 
particle can be obtained by multiplying with (7 — l)/p^ _1 . 
While the specific entropy of individual particles may go 
down if they give up some of their heat energy, the total en- 
tropy of the system increases due to this process, which can 
be interpreted as providing the necessary mixing entropy. 

One important difference of our par amete r izatio n of 'ar- 
tificial heat conduction' to the model of IPricel |2008l) is that 
the mixing only occurs in shear flows, and that contact dis- 
continuities without shear are hence not affected. 

We note that due to the parabolic character of the dif- 
fusion problem, it can be problematic to integrate equation 
(|26|l with an explicit time integration scheme, since the von 
Neumann criterion imposes relatively small timestep limits. 
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where we estimated min 
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L. The CFL criterion is more restrictive than this timestep 
as long as 
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which is usually the case and hence not too restrictive. In- 
deed, so far we have not encountered problems with the 
explicit time integration scheme that is implemented for the 
mixing at present. If needed, an implicit scheme with per- 
fect stabilit y could however be easily adopted, like the one 
discussed in IPetkova fc Springe! (|2009l ). 
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3 ISSUES OF CELL REGULARITY 

A common feature of particle hydrodynamic schemes is 
their ability to automatically provide an adaptive resolu- 
tion. As a result, dense regions are modelled with better 
accuracy thanks to their smaller mean distance of parti- 
cles. But besides the particle number density the regularity 
of the Voronoi cells is an important factor in determining 
the achieved precision, as shown in Appendix IA3I where 
we give quantitative results for the accuracy of our gradi- 
ent estimates as a function of the shape distortions of cells. 
Highly irregular, sliver-like Voronoi cells may also lead to 
very small, computationally costly timesteps, because the 
permissible timestep size is effectively proportional to the 
distance to the nearest neighbour. 

Connected to this problem is the issue of how to safely 
prevent inter-particle penetrations, which is required for a 
proper representation of the fluid with its single valued ve- 
locity field. If two particle approach each other rapidly, it 
is possible that the particles pass through each other unless 
this is prevented with a sufficiently strong artificial viscos- 
ity. If the Voronoi mesh is very irregular and features a large 
number of close particle pairs, it becomes more difficult to 
ensure this, simply because rather large viscous forces that 
act over short timescales are required to prevent the small 
particle separations from becoming still smaller. These prob- 
lems are significantly alleviated if the tessellation is rela- 
tively 'regular', i.e. if cells have a small aspect ratio, and 
if their generating points lie close to the centroids of the 
corresponding cells. 

Figure [2] illustrates another important feature of 
Voronoi meshes, which we may perhaps call 'mesh degen- 
eracy'. In this example, the mesh for two different point 
distributions is shown, but in both cases an identical Carte- 
sian Voronoi mesh results, hence the density and pressure 
estimates are both equal. In fact, one can continue to move 
the groups of four points around the mesh vertices in a mir- 
rored fashion arbitrarily close towards the corners of the 
mesh, without changing the situation. If one now imagines 
adding some random velocity field to the points when they 
are very close, it is clear that it will be much harder to pre- 
vent an erroneous particle crossing in the case where the 
points are far from the centres of their associated cells than 
in the case where they sit right at these centres. 

One might argue that situations as shown in Figure [2] 
are artificial and hence do not affect the simulation of flows 
where Voronoi diagrams are seldom that regular. But situa- 
tions occur where the volume is only slightly increased when 
particles approach. Also considering a finite timestep these 
particles' resistance against a clumping or even interpene- 
tration is then uncomfortably weak. 

If possible, it would therefore be desirable to formulate 
the dynamics such that the mesh automatically maintains a 
certain degree of regularity. We note that this cannot be ex- 
pected to happen by itself for the density estimation scheme 
we implemented thus far. This is again made clear by the 
example in Figure where the pressure gradient vanishes 
in both cases if the specific entropies and masses of all parti- 
cles are the same. Below we therefore consider two possible 
approaches to introduce small correction forces into the dy- 
namics with the goal to rearrange the points to achieve a 
more regular tessellation. 



Figure 2. Two point distributions and their corresponding 
Voronoi tessellations. The important point illustrated by this 
example is that different point distributions may have identical 
Voronoi tessellations. 



3.1 Viscous forces that help to improve order 

Experience with our new VPH scheme shows that especially 
in highly irregular, turbulent flows and in situations with 
strong gravitational forces some cells can become quite ir- 
regular. In this context, we loosely define an irregular cell as 
one whose generating point is substantially displaced from 
the centroid of the cell and/or whose aspect ratio is quite 
high, i.e. a cell with a comparatively large ratio of surface 
area to volume. In this subsection we consider a scheme 
where the artificial viscosity is modified such that it serves 
a second purpose, namely to have the tendency to make cells 
'rounder'. 

To this end we introduce the notion of a 'partial pres- 
sure' for each of the pyramids that make up a Voronoi cell. 
These pyramids are spanned by the Voronoi faces and the 
defining point of the cell, which acts as their apex. We define 
the 'partial pressure' of the pyramid of cell i facing cell j in 
terms of its volume Vy = AijRij/6 (or Vy = in 
2D) and by assigning a share my = rriiAij /^2 k Aik of the 
cell's mass to the pyramid, i.e. its mass fraction is taken to 
be proportional to its contribution to the total surface area 
of the cell. This yields 

The idea is now to define an additional viscous force between 
a pair of particles arising from the difference of this partial 
pressure to the full pressure of the cell. This will lead to rear- 
rangements of the points until the differences in the partial 
pressures of each pyramid to that of the cell become small, 
which happens when the point is approximately equidistant 
to all the faces of its cell, implying a regular cell shape. 
We hence make the ansatz 

(/«**)« = -^(Pij - P + P jt - P 3 )^- (30) 

for 'ordering forces' between a pair of particles, with the 
total force on particle i being given by 

mi(a ovdcr )i = 2j(/ order )ij- (31) 

Here n is a dimensionless parameter describing the strength 
of the effect. These forces are antisymmetric, and in general 
can be both of repulsive and attractive nature. In order to 



8 S. Hefi & V. Springel 



maintain total energy conservation, the work of these forces 
needs to be balanced in the evolution of the entropies of 
the particles, similar to what is done for the ordinary arti- 
ficial viscosity. This yields an additional contribution to the 
entropy change of the form 



is elongated. Specifically, we adopt as Lagrangian 



dsi 
~dt 



(7 - 1) 1 
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Note that a small local decrease of the entropy sum can 
result in principle if order is restored in the particle distri- 
bution, but this effect is small and has played no role in 
all our tests. In refinements of the above scheme, it is also 
possible to make k spatially and temporarily variable. For 
example, we have typically used this scheme in a combina- 
tion with a switch that only sets k > if there is a strong 
local compression, because that is where cell shapes distort 
the most. 

When discussing results, we will refer to this method as 
'partial pressure ordering' or PPO, whereas the method for 
improved cell regularity discussed in the next subsection will 
be referred to as 'shape correction forces'. If none of these 
additional schemes to regularize cell shapes is employed, we 
simply refer to the method as the 'plain Voronoi scheme'. 



3.2 Imposing regularity through the fluid 
Lagrangian 

If irregular cell shapes occur, we ideally would like that small 
adjustment forces appear naturally that tend to make the 
mesh more regular again. These adjustment forces should 
preserve the energy and momentum conservation of the 
scheme. This will automatically be the case if they are de- 
rived from a suitably defined Lagrangian or Hamiltonian. 
We are hence led to modify the fluid Lagrangian slightly 
to include factors that penalize highly distorted cell shapes. 
The idea is that such distorted cells should raise the esti- 
mate of the inner energy slightly, such that they become 
energetically disfavored. 

We consider two ways to measure shape distortions of 
cells, which may either be used individually, or combined. 
One is based on the displacement of a point from the cen- 
troid of its associated cell. The idea here is that it is ad- 
vantageous if a point stays close to the center of a cell. In 
particular, as we will discuss in more detail in Section [5.31 
it turns out that the ordinary VPH scheme is not able to 
support waves in regular Cartesian grids at the Nyquist fre- 
quency, a deficit that could be cured if there is always a 
(weak) restoring force if a point is displaced from the centre 
of its Voronoi cell. 

The other is to measure the shape directly, and to steer 
the particle motion such that high aspect ratios are avoided. 
We construct a shape measure based on the second moment 
of the cell, which we compare to a suitably defined cell ra- 
dius. This measure will have a minimum for 'round cells', 
while severe distortions from roundness (like highly elon- 
gated cells) should trigger restoring forces. 

Both of the above measures of cell regularity can be in- 
troduced into the fluid Lagrangian by multiplying the ther- 
mal energy with correction factors that increase the energy 
slightly if a point is displaced from the centroid, or if a cell 



L = 
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Here is the coordinate of a point, Vk is the volume of its 
corresponding Voronoi cell, and Sk is the cell's centroid. P k is 
the ordinary pressure of the cell, where we set pu — m^/Vk as 
usual. The coefficients /3q and /3i are introduced to measure 
the strength of the correction forces associated with offsets 
from cell centres, or with high aspect ratios, respectively. 
For /3o = /3i = 0, the ordinary fluid Lagrangian of the VPH 
scheme is recovered. 

We define the centroid of a cell as 



Sk 



V k 



rxk(r)dr, 



(34) 



where Xfe is the characteristic function of cell k, i.e. Xk( r ) = 
1 if the point r lies in the cell k, and Xk( r ) = otherwise. 
The shape of a cell is measured via the second moment 



1 

Vk 



(r - s k ) 2 Xk(r) dr. 



(35) 



The factors in squared and curly brackets in the Lagrangian 
of equation (|33p are the adopted energy correction factors 
that raise the energy of distorted cells. Here d counts the 
number of dimensions, i.e. d — 2 for 2D and d = 3 for 
3D. The factor Vf is hence proportional to the 'radius' 



vV d of a cell squared. /?o measures the strength of 
the effect of displacements of points from the centroid of a 
cell, while /3i is the corresponding factor for the aspect-ratio 
factor. The constant fiz is only introduced to prevent that 
even round cells lead to a significant enhancement of the 
thermal energy. For perfectly round cells, we expect in 2D 
approximately circles for which w 2 . = V 2 ^ d /(2n), hence we 
pick P2 — l/(27r). In 3D, we have spherical shapes instead 
and we pick /3 2 = 3/5(3/4vr) 2/3 . 

The equations of motion for the Lagrangian (|33l) can be 
derived in closed form for the Voronoi mesh, but due to the 
length of the resulting expressions we give their derivation 
in Appendix [B] The advantage of using the Lagrangian to 
obtain the cell-shaping forces is that the scheme then still 
accurately conserves total energy, momentum and entropy, 
while at the same time remaining translational and rotation- 
ally invariant. Also, the correction forces are 'just right' to 
achieve the desired regularity of the mesh, something that is 
difficult to achieve with any heuristic scheme to derive such 
forces, like the one we tried in the previous subsection. 

In our results section, we show that this method is in- 
deed capable of maintaining nicely regular meshes in the 
sense described above. However, we also caution that the 
stronger the extra forces are, the more unwanted features 
start to appear as well. First of all, the extra forces may 
introduce subtle deviations from the dispersion relation of 
an ideal gas, and may lead to spurious motions in situations 
with pressure equilibrium. The second, probably more se- 
rious side effect of this method may occur when the cells 
cannot easily relax to the desired regular cell structure, for 
example along a strong jump in density. In this case a pres- 
sure anomaly may develop due to the cell-shaping forces, 
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similar to what is found in SPH across contact discontinu- 
ities. Still, we find that moderate values of /3o and Pi help to 
improve the accuracy of the VPH scheme without distorting 
the inviscid dynamics of an ideal gas too much. 



4 IMPLEMENTATION 

We have implemented the above hydrodynamical parti- 
cle model into the cosmological TreeSPH simulat ion code 
GADGET-3, an updated v ersion of GADGET-2 (|Springell 
2005; ISpringel et al. I 120011 ). This code is parallelized for 
distributed memory machines, and offers high-performance 
solvers for self-gravity as well as individual and adaptive 
timestepping for all particles. Our strategy in our modifica- 
tions has been to implement Voronoi-based particle hydro- 
dynamics as an alternative to SPH within the GADGET-3 
code. This is, in particular, ideal for facilitating comparisons 
between SPH and our Voronoi-based scheme, and it also al- 
lows us to readily use all the non-standard physics already 
implemented in GADGET-3 (e.g. radiative cooling and star 
formation) for calculations with Voronoi-based fluid particle 
dynamics. 

The primary new code needed in GADGET-3 is an effi- 
cient mesh construction algorithm. To this end we adapted 
and modified the par allel Delaunay triangulation engine 
from the AREPO code l|Springel|l2009l ). and turned it into an 
optional module of the SPH code GADGET-3. In brief, the 
tessellation code uses an incremental construction algorithm 
for creating the Delaunay tessellation. Particles are inserted 
in turn into an already existing, valid tessellation. To this 
end, in a first step the tetrahedron in which the new point 
falls is located, and then it is split into several new tetrahe- 
dra, such that the inserted point becomes part of the tetrahe- 
dralization. However, some of the new tetrahedra may then 
not fulfill the empty-circumsphere property, i.e. the tessel- 
lation is not a Delaunay triangulation any more. Delaunay- 
hood is restored in a second step by local flip operations 
that replace two adjacent tetrahedra with three tetrahedra, 
or vice versa, until all tetrahedra fulfill again the empty- 
circumsphere property. At this point, the next particle can 
be inserted. We have implemented the mesh construction 
both for 3D and 2D within the GADGET-3 code and paral- 
lelized it for distributed memory machines. 

At the beginning, a large tetrahedron is constructed 
that encloses the full computational domain. The boundary 
conditions (always adopted as periodic at the moment) of 
the rectangular computational domain as well as the bound- 
aries arising from the domain decomposition are treated 
with 'ghost particles'. One technically difficult aspect is to 
make the tessellation code completely robust even in the case 
of the existence of degenerate particle distributions, where 
more than 4 points lie on a common circumsphere (or more 
than 3 points are on a common circumcircle) . Detecting such 
a case robustly and correctly in light of the finite precision 
of floating point arithmetic is a non-trivial problem. How- 
ever, the incremental insertion algorithm requires consistent 
and correct evaluations of all geometric predicates, other- 
wise it will typically fail in situations with degeneracies or 
near-degeneracies. We solve this problem by monitoring the 
floating point round-off in geometric tests, and by resorting 



to exact arithmetic in case there is a risk that the result of 
a geometric test may be modified by round-off error. 



5 TEST RESULTS 

In this section, we discuss a number of test problems carried 
out with our new hydrodynamical particle method, focus- 
ing in particular on regimes where differences with respect 
to SPH can be expected. An application of the method in 
full cosmological simulations of galaxy formation will be pre- 
sented in future work. 

5.1 Surface tension 

In standard entropy-conserving SPH with particles, there is 
a subtle surface tension effect across contact discontinuities 
with a large jump in density. This can be understood as a 
result of the desire of SPH to suppress mixing of the two 
phases, because this is energetically unfavorable for fixed 
particle entropies. For the mixed state, approximately the 
same average density would be estimated for each particle, 
which leads to a higher estimate of the thermal energy, un- 
less the thermodynamic entropies are averaged between the 
particles as well, which is an irreversible process in which en- 
tropy is in fact produced if the total energy stays constant. 

To demonstrate the existence of the surface tension ef- 
fect, we have prepared (in 2D for simplicity) an overdense 
ellipse in a thin background medium, at pressure equilib- 
rium. For definiteness, the density of the ellipse was set to 
P2 — 4, that of the background medium to pi = 1, with 
a pressure of P = 2.5. 3854 equal mass particles in a peri- 
odic box of unit length on a side were used to set up the 
experiment. The particles have been arranged on a coarse 
Cartesian grid in which an ellipsoidal region was excised. 
This region has then be filled with a finer Cartesian grid. 
The specific entropies of the two sets of particles were ini- 
tialized differently such that an equal pressure for the two 
phases results. No attempt was made to somehow soften the 
transition between the two phases. Figure [3] shows the ini- 
tial configuration, as well as the particle distribution after 
a time t — 7, both for SPH and for the Voronoi-based fluid 
particle approach. 

Even though the pressures of the particles are formally 
equal for all particles in the initial conditions, the ellipse 
transforms to a circle when SPH is used. In contrast, for 
the VPH scheme, the same experiment maintains the initial 
shape of the ellipse, modulo some small rearrangements of 
the points near the boundary, since the initial set-up was not 
in perfect equilibrium (due to the fact that the point distri- 
butions of the two Cartesian grids used to set up the two 
phases do not match seamlessly at the boundary). Clearly, 
the numerical realization of the contact discontinuity in SPH 
gives rise to a spurious surface tension, and this in turn will 
suppress Kelvin-Helmholtz instabilities below a certain criti- 
cal wavelength. The VPH approach does not have this prob- 
lem and can in principle accurately support a contact dis- 
continuity at each face boundary between individual cells. 
It needs to be stressed however that also in the Voronoi 
scheme no mixing of the entropies at the particle level hap- 
pens. If the particles of two phases were simply spatially 
mixed while keeping their specific entropies constant, the 
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resulting medium would not be at a single temperature or 
density. We note that the effect is also present for 
a set-up with unequal particle masses. In this case, 
the variationally derived entropy-conserving version 
of SPH will however lead to a change of the effective 
number of neighbours across the contact discontinu- 
ity, because it keeps the mass in the kernel volume 
constant. 

5.2 Sod shock tube 

The classic Sod shock tube tests examine the ability of a 
hydrodynamic scheme to reproduce the basic wave struc- 
ture that appears in the Riemann problem, namely shock 
waves, contact discontinuities and rarefaction waves. Also, 
comparison to the analytic solution gives a useful quantita- 
tive benchmark for the accuracy of a scheme. 

We consider gas that is initially at rest. In the left half- 
space, the pressure is Pi = 1.0 and the density is pi — 1.0, 
whereas in the right half-space we adopt P2 = 0.1795 and 
P2 — 0.25. The adiabatic index is set to 7 = 1.4. The same 
sod shock parame ters have previously been used in a number 
of code tests (e.g. Hernquist fc Katz 19891 ; iRasio fc Shapiro! 
Il99ll ; IWadslev et al.ll2004ISpringe]||2005l ). When the evolu- 
tion begins, a shock wave of Mach number M — 1.48 travels 
into the low-pressure region, and a rarefaction fan moves 
into the high pressure region. In between, a moving contact 
discontinuity develops. 

In our numerical test of this problem with the VPH 
scheme we use a 3D setup in a box with dimensions (20, 1, 1), 
where the left and right halves are filled with particles ar- 
ranged on a Cartesian grid. Altogether 8370 particles with 
equal masses were used as initial condition. The evolution 
was then carried out with our default settings for the artifi- 
cial viscosity until t = 3. In Fig. [4j we compare the numerical 
result to the analytical solution at this time. Reassuringly, 
we find quite good agreement of the VPH scheme with the 
analytical solution. In comparison to SPH, the rarefaction 
wave in the Voronoi simulation shows a slight dip at the low 
density end, but not the hump at the high density end that 
is typical of SPH. Another difference is that the contact 
discontinuity is sharper and better preserved in the VPH 
approach. All of the described features appear to be largely 
independent of how the points are distributed initially as 
long as the gas is relaxed on both sides. In particular, the 
arrangement on a Cartesian grid does not lead to any notice- 
able artifacts compared with simulations where the initial 
point distribution is less ordered and has a glass-like config- 
uration. However the results are somewhat less accurate for 
irregular grids (see also lA3|) when the plain Voronoi scheme 
is used. 

We have also examined the influence of the extra forces 
that can be enabled to improve the regularity of the mesh 
(see Sections l3.ll and l3.2p . To this end we tested the effect of 
the PPO scheme as well as the shape correction method with 
parameters up to /3o = 1-2, /3i = 0.1. When these ordering 
forces are invoked, the results for the sod-shock test are in 
general not influenced much, but the particle noise around 
the analytical solution is reduced. We also found that the 
extent of postshock oscillations for weaker artificial viscosity 
settings tends to be reduced if these ordering methods are 
used. 
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Figure 4. 3D Sod shock tube simulation at t = 3.0. We compare 
results from our Voronoi particle scheme (blue crosses) and the 
SPH result with 32 neighbours (red diamonds) with the analytical 
solution (solid lines) in terms of the density, pressure, entropy 
profiles and internal energy. 



5.3 Dispersion relations 

Even though this may seem like a simple test, it is actually 
important to check how well our new method can simulate 
small-amplitudeQ acoustic waves, especially at low resolu- 
tion when few points per wavelength are available. We are 
especially interested in how accurately the expected disper- 
sion relation is reproduced in this regime, i.e. whether such 
waves propagate with the correct speed of sound. A sec- 
ondary question is how strongly such waves are damped by 
the artificial viscosity in the scheme. 

To measure the dispersion relation, we set up small- 
amplitude standing waves in a periodic box and measure 
their oscillation frequency. In Figure [S] we compare re- 
sults for SPH with our Voronoi-based fluid particle model, 
with and without shape correction forces, as a function of 
wavenumber. The wavenumber is normalized to the Nyquist 
frequency of the initial particle grid, such that fc/^Nyquist = 1 
corresponds to the shortest wave that can be represented by 
the particles. In this standing wave, neighbouring particles 
oscillate 180 degrees out of phase 'opposite' to each other. 

A first important result made clear by Figure [5] is that 
for the standard VPH scheme the oscillation frequency for 



1 The amplitude needs to be small in order to prevent wave steep- 
ening. 
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Figure 3. Surface tension effect in SPH. The left column shows an overdense ellipsoidal region shortly after it is set-up at t = in 
pressure equilibrium within a thinner background. When evolved with SPH, the ellipsoid slowly transforms into a sphere, as shown by 
the state of the system after time t = 7 (top left). In contrast, the Voronoi scheme can preserve the shape much better (bottom panels) 
and shows no sign of surface tension effects. 



fc/^Nyquist = 1 drops to zero, or in other words, such waves 
are not supported by the scheme at all. This is readily un- 
derstood from the degeneracy effect pointed out in Figure [2] 
If particles are set-up such that they 'collide' in a pairwise 
fashion, then there is nothing in the reversible part of the 
dynamics of the VPH scheme that can prevent an interpar- 
ticle penetration, simply because the pressure gradient stays 
zero in this case. However, this situation is exactly the one 
encountered if we prepare a standing wave at the Nyquist 
frequency of an initially regular particle grid. The wave will 
not oscillate since the pressure gradient will remain zero, 
and therefore particle crossings would be inevitable (unless 
prevented by the artificial viscosity). This is potentially a se- 
rious shortcoming of the VPH scheme in its standard form, 
as it means that it cannot treat waves at around the Nyquist 
frequency properly. 

However, the shape correction force due to j3o has ex- 
actly the right property to make these small waves oscillate 



again. In fact, we can calculate what value of fio is required 
to reproduce the dispersion relation at fc/fcNyquist = 1 ex- 
actly. For this value of fio ~ 5.5, we however also get slightly 
too stiff behaviour of the fluid for somewhat longer wave- 
lengths, as shown by Fig. [5] A value of around /3o ~ 3 rep- 
resents a good compromise, and in particular yields a more 
accurate dispersion relation than SPH for all k. We also note 
that for certain numbers of neighbours, the SPH result is in- 
accurate at all wavelength; here the numerical soundspeed 
shows an offset relative to the expected sound speed, which 
is presumably a result of a bias in the density estimate for 
the background density. 



5.4 Density noise and regularity in a settled 
particle distribution 

In this subsection we want to examine the level of noise 
present in a relaxed region of gas of constant specific entropy, 
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Figure 5. Dependence of the numerical sound speed on 
wavcnumber, expressed in units of the Nyquist frequency of the 
underlying particle grid. For the standard VPH scheme, waves 
at the Nyquist frequency are not propagated properly (hence no 
frequency can be measured at k = k^S), but this is remedied 
by the shape correction forces. If they are invoked, the resulting 
dispersion relation becomes more accurate than that of SPH for 
all k. 

as it may arise somewhere within a larger, self-consistent 
simulation. To mimic this situation, we start from a distribu- 
tion of points arranged on a Cartesian grid and impose a ran- 
dom Gaussian velocity field with dispersion (v 2 ) — 0.05 (? s 
and zero mean. The idea is that these velocity fluctuations 
break the initial grid symmetry and will then get damped 
away by the artificial viscosity, which is here set to a high 
value to speed up the process of settling to a new pressure 
equilibrium. Since we want to retain the initially equal val- 
ues of the specific entropies per particle, we disable the en- 
tropy source term for the viscosity in this experiment. Once 
the new equilibrium for an irregular particle distribution is 
achieved, we can then examine the noise properties of this 
particle representation of a constant density, constant pres- 
sure gas. 

For SPH with N = 16 neighbours, we find that the 
particles settle into several domains in which the points are 
quite regularly distributed, based on visual inspection. The 
estimated density values pi for the particles are not all equal 
though, instead they show a distribution with rms-scatter 
equal to ~ 1.4%, and also a small bias relative to the ex- 
pected value equal to the mean density (p) — Nm/V of the 
full volume. 

In contrast, the standard VPH approach creates a dis- 
tribution in which the density values are essentially single- 
valued, and are all very close to (p). This means that the 
cells have all equal volume, and the residual pressure fluctu- 
ations, if any, are extremely small. However, the geometry 
of the Voronoi tessellation is quite irregular and features 
numerous cells with relatively large aspect ratios, or with 




Figure 6. Final mesh geometry in VPH in a 2D settling test, 
carried out without (top) or with shape correction forces (bottom) 
based on O = 1.0 and ft = 0.01. 



points close to cell boundaries. This can be seen in Figure [6] 
where we show a plot of the final mesh for a test case carried 
out in 2D. 

It is now interesting to repeat the test for the case when 
shape correction forces according to Section|3]2]are included. 
As desired, the final mesh becomes much more regular in 
this case, as seen in the corresponding example included 
in Figure [6] However, even in the final equilibrium state 
the correction forces do not necessarily completely vanish 
in this case. Instead, they are compensated by small resid- 
ual pressure (and hence also density) fluctuations. This is 
demonstrated in the distribution functions of the density 
values for the three cases we considered, which we give in 
Figure [7] Here the shape correction case yields a somewhat 
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Figure 7. Density distribution functions of the particles in a 
settling test for constant entropy gas, carried out with different 
schemes. We compare SPH (red), ordinary VPH (black), the PPO 
version of VPH (blue), and VPH with shape correction forces 
based on j3o = 1.0 and /3i = 0.01 (in green). The distributions 
were measured at a time when the initial kinetic energy had de- 
cayed to E kin ft 0.001B kin (t = 0). 




Figure 8. Cell-regularity of a noisy flow after relaxation. To char- 
acterize the regularity of the cells, we simply consider the distri- 
bution of the distance of the points to their cell's centroids, in 
units of the 2D cell radius r = y^V/n. The black histogram shows 
the distribution for the ordinary Voronoi scheme, blue shows the 
PPO version, and green lines give the result for the Voronoi with 
shape correction forces derived from the Lagrangian. The distri- 
butions were measured at a time when the initial kinetic energy 
had decayed to _B kin ft 0.001E kin (t = 0). 



broader distribution, similar to SPH, but without showing 
a bias. 

In Figure [8] we give a quantitative measure for the cell- 
regularity (here taken as the distribution function of the 
normalized displacement of points from the centres of their 
cells) of the final meshes in the Voronoi-based simulations. 
Even moderate values of the coefficients /?o and /3i can dras- 
tically improve the regularity of the particle distribution 
while introducing less noise in the density estimate than 
anyway present in SPH. 

5.5 Point explosion 

If energy is injected at a point into a cold gas at constant 
density, a spherical blast wave will develop. The Taylor- 
Sedov solution provides an analytic solution for this self- 
similar problem, which is a useful test involving very strong 
shocks. We have set-up this problem in 3D, using unit 
background density (represented with a Cartesian mesh), 
7 = 5/3 and vanishingly small initial specific entropy com- 
pared to the injected energy of E = 1. We inject the energy 
into the centre of the domain at time t = 0. To avoid that 
the evolution is strongly affected by the non-spherical ge- 
ometry of the central Voronoi cell, we have spread out the 
energy with a Gaussian kernel with a radius of about 4 mean 
particle spacings. 

We note that this set-up can be especially sensitive to 
the problem of particle crossing when a too low viscosity 
and individual timesteps are used. In the latter case, the lo- 
cal Courant timestep of particles outside of the explosion is 
initially very big. When the supersonic shock front arrives, 
such particles may then still live on a too large timestep, 



such that they are effectively overtaken by the shock, creat- 
ing severe artifacts in the evolution of the shock front. We 
have addressed this in our test by imposing a low enough 
maximum timestep for all particles, but more sophisticated 
schemes to set the timesteps, which guarantee that it is re- 
duced before the shock arrives, can of course be impleme nted 
in principal (see lSaitoh fc Makinoll2(5o9l ; ISpringelll2009D . 

In Figure we show the radial density profile at t — 
0.04 for different resolutions corresponding to 2 x 32, 
2 x 64 and 2 x 128 particles, and compare to the expected 
analytic solution. The expected solution is captured reason- 
ably well, with a similar quality as in SPH codes, which 
is shown for comparison in the low-resolution case. 
The VPH result converges nicely to the analytic so- 
lution as the resolution is increased. In particular, the 
shock location is well reproduced, albeit with a small pre- 
shock increase of the density. When shape correction forces 
are added as described in Section \3. 21 only negligible differ- 
ences in the overall quality of the result are found, but the 
scatter around the azimuthally averaged solution is reduced. 
In Figure 11 L we check the energy conservation in 
the blast wave problem, both for the VPH and SPH 
simulations at the 32 3 resolution. The energy error 
is negligibly small, as desired. The high-frequency 
oscillations in the total energy of the SPH run stem 
from the small but finite jumps of the SPH smooth- 
ing lengths from timestep to timestep, an effect that 
is absent in the VPH simulation. 
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Figure 9. Sedov- Taylor point explosion problem, calculated in 
3D with the basic VPH scheme. From top to bottom, 
we compare the radial density profile at t = 0.04 with 
the analytical solution (solid line), calculated with 32 3 , 
64 3 , or 128 3 particles, as labeled. In the top panel, we also 
compare the VPH results (circles) with SPH (diamonds). 
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Figure 10. Time evolution of the thermal and kinetic energies in 
the Sedov- Taylor point explosion problem, simulated with SPH 
(black dashed line) and the Voronoi scheme with shape correc- 
tion forces (red solid line). The bottom panel compares the total 
energy error in the two schemes. 



5.6 Kelvin-Helmholtz instabilities 

Kelvin-Helmholtz (KH) instabilities occur in regions of 
strong shear, which is especially common at contact discon- 
tinuities between two fluid phases. An initially small trans- 
verse perturbation along the interface becomes amplified 
and grows in linear theory according to oc exp(i/iicH)- After 
a few characteristic timescales 



pi + P2 
2kv^/p 1 p 2 ' 



(36) 



an initially wave-like perturbation becomes large and non- 
linear, developing the typical KH-rolls. Here p\ and p2 are 
the densities of the two media, v is the velocity jump par- 
allel to their common interface, and k = 27r/A por t is the 
wavenumber of the perturbation with wavelength A pcrt . For 
an ideal gas, all wavelengths are unstable, and the smallest 
wavelengths grow fastest. 

The KH instability is especially important for the de- 
velopment of turbulence, and is thought to play a promi- 
nent role in stripping and mixing processes occurring dur- 
ing galaxy formation. Recently, a number of studies have 
pointed out that standard SPH has problems to correctly 
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capture the KH instability when the initial conditions con - 
tain sharp density gradients (jAgertz et al.ll2007j ; |Pricef 2008). 
In certain cases, the instability is suppressed completely 
and does simply not grow. This can in part be under- 
stood in terms of the surface tension effect present in 
SPH, as described earlier, because surface tension sup- 
presses the growth of KH instabi lities below a critical wave- 
length (jLandau fc Lifshit d Il966h . Furthermore, the asym- 
metric particle density at the interface causes a rearrange- 
ment of the points in SPH, such that a 'gap' in the sampling 
appears that causes r elatively large error s in the pressure 
forces at the interface l|Agertz et al.ll2007l ). 

It is therefore very interesting to test how well the VPH 
approach does in this respect. Since VPH does not exhibit a 
surface tension effect, it offers the prospect of a better treat- 
ment of the KH instability. In Figure QT] we show results for 
a KH test calculation, carried out with different particle- 
based hydrodynamic schemes. Our two-dimensional initial 
conditions consist of 7 = 5/3 gas with density pi — 2 in 
the stripe \y — 0.5| < 0.25, moving to the right with veloc- 
ity vi — 0.5, and of gas with density p 2 = 1 and velocity 
V2 = —0.5 in the region \y — 0.5| ^ 0.25. The pressure was 
initialized everywhere to P — 2.5, and a periodic domain of 
unit length on a side was used. In total 261760 points sample 
the gas distribution, with a mean spacing of 0.0023 for the 
low-density gas, and 0.0017 for the high density gas. Hence 
the two phases were represented with approximately equal 
mass particles. Note that the initial discontinuity was im- 
posed as a perfectly sharp jump in these initial conditions, 
following previous studies of this problem. We remark how- 
ever that it is somewhat questionable whether such sharp 
jumps are not introducing an inconsistency with the ba- 
sic premises of SPH calculations, which can only represent 
smoothed density fields. In order to seed an initial perturba- 
tion, we imposed a vertical perturbation on the y-positions 
of the form 

Sy(x) = ao sin(47ra;/L), (37) 

where L is the boxsize, and ao = 0.006 is the amplitude of 
the initial perturbation. 

The three columns of Figure [11] compare the results for 
the ordinary VPH scheme (left), the VPH method with the 
additional ordering viscosity of the PPO scheme (middle), 
and SPH (right), at time t = 1.2. From top to bottom, 
we show specific entropy maps, pressure maps, the particle 
distribution, and vorticity maps. All maps were here gen- 
erated by linearly interpolating a Delaunay tessellation of 
the points, allowing to extend the points' properties as read 
from the simulation files to continuous fields. In addition to 
the maps, we show in Figure [12] the growth of the instabil- 
ity until time t = 1.2, quantified in terms of the amplitude 
of the seeded velocity mode as measured in the Fourier- 
transformed Vy-Qeld. 

We see right away that the VPH scheme captures the 
KH instability best. Its primary KH billow has evolved fur- 
thest, and it triggered the growth of smaller-scale secondary 
billows. In contrast, the SPH result shows only an anemic 
growth of the instability. In the SPH pressure map a strong 
pressure anomaly is visible at the interface. This surface 
effect effectively suppresses all small-scale KH-instabilities, 
and the two phases stay separate because the associated sur- 
face tension suppresses a breaking up of the interface. As a 



result, the instability cannot cascade down to smaller scales. 
The PPO scheme shown in the middle column lies literally 
in the middle in this respect. The growth of the primary 
KH mode is very similar to that found in the VPH scheme. 
However, the additional viscosity introduced in this scheme 
to produce highly regular cells substantially attenuates the 
growth of secondary small-scale KH instabilities. The same 
effect can be seen Figure [12] for the case with shape correc- 
tion terms, whereas additional heat diffusion (according to 
Section 12. 5[) does not affect the growth rate of the excited 
mode in this simulation. We note however that the smaller 
growth depends on the strength of the additional ordering 
forces that are invoked. The result shown here was calcu- 
lated with k = 1, which we consider the maximum that one 
may ever want to use. For more reasonable smaller values, 
intermediate results that are close to that of the plain VPH 
scheme are obtained. 

Interestingly, the vorticity V x v in the ordinary VPH 
scheme is clearly largest overall, especially for the larger 
modes and in the central region of the primary KH billow. 
Here the rotation is so fast that the pressure shows a no- 
ticeable depression, which counteracts the centrifugal forces 
from the rotation with pressure gradients. Unfortunately, 
the ability of the pure Voronoi simulation to sustain vor- 
tices comes at the expense of a larger particle irregularity. 
As the enlargement with the particle distribution shows, the 
particles tend to form Voronoi cells with quite high aspect 
ratios, reducing the accuracy of the gradient estimates and 
requiring relatively large settings for the artificial viscosity 
to prevent interparticle penetrations. In contrast, the PPO 
variant of the Voronoi scheme produces highly regular parti- 
cle spacings, and in this respect resembles SPH. However, in 
this example calculation with k = 1, the scheme then devel- 
ops effectively a much higher intrinsic shear viscosity, which 
tends to transform the differential rotation into a rotation 
of numerous ordered domains. The shear of these domains 
shows up in the differential vorticity maps in a distribute 
fashion. Once the damping of the differential rotation gets 
too strong, the primary vortex is affected as well, as seen in 
the reduction of the central pressure gradient in the PPO 
scheme. 

Finally, we test the shape correction forces discussed in 
Section T3. 21 and our scheme for the treatment of subresolu- 
tion mixing introduced in Section [2. 51 with the same initial 
conditions. We show in Fig.[T3]the resulting entropy maps of 
two further calculations of the KH instability test. In both 
panels, we show runs of the VPH scheme that make use 
of the shape correction forces derived from the Lagrangian, 
with f3o = 1.2 and /3i — 0.1. In the bottom panel, we have in 
addition activated the artificial heat conduction due to local 
shear with \ = 0.25, which models mixing of the fluids at 
the scale of the resolution. The two variants are qualitatively 
similar, but the artificial heat conduction has clearly washed 
out the sharp discontinuity in the entropy at the fluid inter- 
face. This resembles more closely the results of mesh-based 
finite-volume hydro codes. Compared to the results of pure 
VPH in Figure [TT] we see that the scheme with shape cor- 
rection forces has also an effectively enlarged viscosity, quite 
similar to the PPO approach. Overall, it is clear that even 
without a subresolution mixing model the Voronoi based 
particle hydrodynamics does significantly better in the KH 
test than standard SPH. 






Figure 11. Simulations of the KH-instability with different particle-based methods. The left column shows the plain Voronoi scheme, 
the middle column Voronoi with PPO, and the right column SPH. From top to bottom, we show maps of specific entropy, maps of the 
pressure, the point distribution, and maps of the differential rotation V X v. The maps show x = [0.f8, 0.58], y = [0, 0.4] of the periodic 
simulation domain, while for clarity the particle distribution is shown only for x = [0.28,0.48],?/ = [0.f,0.3]. 
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Figure 12. Growth rate of the KH-instability for the same initial 
conditions as in Figure [TT] We show the amplitude of the seeded 
mode in the velocity field, measured by Fourier-transforming the 
v y field. We give results for standard VPH (green), for the Voronoi 
scheme with additional shape correction forces and the prescrip- 
tion for mixing discussed in Section 12.51 (dotted blue) , and for 
SPH (red). The dashed black lines indicate the slope expected 
for an exponential growth of the instability according to linear 
theory 




5.7 The 'blob test': mass loss of a gas cloud in a 
supersonic wind 

A challenging test problem for hydro dynamical codes has 
been proposed by lAgertz et al.l |2007l ) . The setup consists 
of an overdense spherical cloud in pressure equilibrium with 
the surrounding hot medium. This background gas is given a 
large velocity, so that the cloud feels it as a supersonic head 
wind. The test is motivated by astrophysical situations such 
as the stripping of gas out of the halos of galaxies as they fall 
into larger systems. It is a three-dimensional problem that 
involves many different non-linear hydrodynamical phenom- 
ena, including shocks, Kelvin-Helmholtz instabilities, mix- 
ing, and the generation of turbulence. Because of this com- 
plexity, an analytical solution for the problem is not known. 
The general expectation is that the wind will compress the 
cloud, accelerate it, and strip some of its gas by developing 
fluid instabilities as it streams past the cloud. 

In terestingly, in the test calculations of lAgertz et al.l 
(2007), substantial differences were found in the mass loss 
rates of the cloud when calculated with Eulerian mesh codes 
and with SPH. Whereas the mesh codes led to an eventual 
complete destruction of the cloud, the mass loss rate was in 
general smaller in SPH, such that some cloud material still 
remained once the partially destroyed cloud was accelerated 
to the wind speed and the mass loss stopped. Given these 
qualitatively different outcomes, it is interesting to test how 
the VPH scheme performs on this problem. 

We used the s ame initial conditions as employed by 
lAgertz et al.l ((2007), in the version with 10 6 particles. The 
setup consists of a periodic box with extension [0, 2000] x 
[0, 2000] x [0, 8000] kpc. The background 'wind' gas has den- 
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Figure 13. KH-instability test simulated with the Voronoi 
scheme with additional shape correction forces, based on the same 
initial conditions as in Figure ITT1 and again at time t = 1.2. The 
maps show the entropy distribution without (top) and with (bot- 
tom) additional heat diffusion terms to model subrcsolution mix- 
ing, as described in Section 12.51 

sity Pwind = 4.74 x 10~ 34 gcm -3 , temperature T w ind = 
10 7 K, and a velocity v wind = (0, 0, 1000) kms" 1 . The cloud 
has a radius -Rdoud = 197 kpc, and is placed initially at 
r c ioud = (1000, 1000, 1000) kpc. Its density is 10 times higher 
than the background, p c i ou d = 10p w i n d, while at the same 
time being 10 times colder, T c i ou d = 10 6 K. This yields a 
sound speed of Cwind = 371 kms" 1 for the wind, and an 
expected characteristic timescale of order tkh — 2 Gyr for 
the development of large KH instabilities in the shear-flow 
around the cloud. 

In Figure 1141 we show the time evolution of density 
slices through the central plane of the simulation box, calcu- 
lated with our new Voronoi scheme. In agreement with both 
grid-based and SPH codes, a bow shock is produced ahead of 
the cloud. The cloud gets compressed and accelerated under 
the ram pressure of the wind, and the wind that streams past 
the deformed and slowly accelerating cloud induces Kelvin- 
Helmholtz instabilities at its surface which strip material 
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and produce a turbulent wake. At the stagnation point of 
the flow in front of the cloud, the pressure eventually breaks 
through along the axis of symmetry. The resulting "smoke 
ring" is then still exposed to Kelvin-Helmholtz induced tur- 
bulence while it is being accelerated to the velocity of the 
background stream. Due to the periodic boundary condi- 
tions we note that the bow shock extends past the domain 
boundary and then back inwards again from the other side, 
reaching the cloud at about t « 3 tkh- This leads to a recom- 
pression of the remainder of the cloud which can temporarily 
raise the number of particles that are still counted as cloud 
members. 

The mass loss as a function of time is displayed in Fig- 
ure !15l for different particle-ba sed hydrodynamical schemes. 
We follow lAgertz et"aD (|2007h and consider a particle to be 
still part of the cloud when its density is still larger than 
p > 0.64p c i oud , and its temperature fulfills T < 0.9 T w i n d- We 
show results for four different calculations in total. The stan- 
dard VPH scheme is shown in blue. Interestingly, it leads to a 
complete destruction of the cloud at time t « 3 tkh , a result 
which is actually surprisingly clos e to the high-resolu tion 
mesh-based calculations reported in lAgertz et al](|2007l ). On 
the other hand, the two SPH-based results (shown in red 
and black) calculated with the GADGET2 code do not re- 
sult in a destruction of the cloud. Instead at time t = 5 tkh, 
still about half the mass of the original cloud can be char- 
acterized as residual cloud mat erial. This is even slightly 
larger than what was reported bv lAgertz et alj (|2007T ) for the 
GASOLINE code. We have however found that the formula- 
tion and the strength of the artificial viscosity can influence 
this result significantly. Also, we confirmed that integra- 
tion of the entropy as independent thermodynamic variable 
(which is the default in GADGET2) results in less stripped 
material than when the thermal energy is integrated as done 
in GASOLINE. This is however probably largely a result of 
the particular initial c onditions used here ; the contact dis- 
continuity in the ICs of lAgertz et~ai1 l|200it ). that we employ 
here, has been relaxed using the traditional SPH formulation 
of GASOLINE, creating a pressure blip. When this is then 
used to initialize the entropies integrated in GADGET2, a 
spurious entropy blip is created that further amplifies the 
initial sampling 'gap'. 

When strong shape correction forces are introduced into 
the VPH formalism, we find an intermediate result between 
plain VPH and SPH. In this case, a small 10% remnant of 
the cloud remains at time t = 5 tkh • This is consistent with 
our earlier findings for the KH instabilities. While standard 
SPH can be expected to suppress the KH-instabilities signif- 
icantly at the cloud interface, it tends to underestimate the 
rate of stripping. Our new VPH method does not show this 
problem, but if viscous forces are introduced that guaran- 
tee very regular mean particle separations, some small-scale 
suppression of fluid instabilities can be reintroduced. 

We finally note that the use of a time variable viscosity 
as presently implemented in our code has not changed the 
mass-loss curves significantly. The reason is that the rele- 
vant particle viscosities are pushed to a large value as they 
pass through the bow shock, and stay at large values in 
the complicated flow around the cloud surface. Only at late 
times the mass loss tends to become faster as a result of the 
effectively lower viscosity. 
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Figure 14. Time evolution of the density for a gas cloud in a 
supersonic wind. From top to bottom, we show density maps 
normalized to the initial wind density at times t = 0.75 tkh, 
t = 1.5 tkh i an d t = 2.25tkh i n the central plane of the simu- 
lation box. Here the standard Voronoi scheme with 10 6 particles 
was used. 
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Figure 15. Mass loss of a gas cloud in a supersonic wind (the 
'blob test'), simulated with 10 6 particles and different hydrody- 
namical schemes. We show the remaining mass associated with 
the cloud (particles that fulfill p > 0.64p c i out j and T < 0.9T w ; n( j) 
as a function of time (in units of tkh) f° r different hydrody- 
namical schemes. The different colors refer to the SPH codes 
GADGET2 (black), modified GADGET2 with an GASOLINE- 
like integration scheme (red), standard VPH (green), VPH with 
dynamic viscosity (see Equation 1221 (light blue) and VPH with 
heat diffusion according to Section 12.51 (dark blue) 



5.8 Gravitational collapse of a gas sphere 

Finally, we consider a three-dimensional problem with self- 
gravity, the so-called 'Evrard-collapse' (jEvrardl Il988l ). It 
consists of an initially cold gas cloud, with a spherically 
symmetric density profile of p(r) oc 1/r. The total mass, 
outer radius and gravitational constant are all set to unity, 
M = R = G = 1, and the initial thermal energy per 
unit mass is set to u = 0.05. In this configuration the 
sphere is significantly underpressurized. It hence collapses 
essentially in free-fall, until it bounces back at the centre, 
with a strong shock running through the infalling material. 
The sphere then settles into a new virial equilibrium. As 
this problem involves large conversions of potential gravita- 
tional energy into kinetic energy and thermal energy (and 
back), as well as strong shocks, it is a challenging and use- 
ful test for hydrodynamic codes that are applied to struc- 
ture formation problems. For this reason, it has been widely 
used as a test for a number of SPH codes (e.g. Evrardl 
19881; iHernauist fc Kat3 Il989l; [steinmetz fc Muellerl I l995| 
Dave et al.lll997l ; ISpringel et al.ll200ll ; IvVadslev et all 12004 ; 
Springelll2005r i. 

For our test we create a realization of the sphere by 
stretching a Cartesian grid appropriately, such that the de- 
sired initial density profile is obtained. Because the Voronoi 
scheme needs to tessellate a well-defined total volume, we 
cannot impose vacuum boundaries in the same way as in 
SPH. Instead, we embed the sphere in a box and use a back- 
ground grid of particles with extremely rapidly falling 
density profile outside of the sphere, such that the total 
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Figure 16. Cell regularity without (bottom) and with cell regu- 
larization forces (top) in a simulation of a sphere of gas collaps- 
ing under self-gravity (Evrard collapse problem). In the bottom 
panel, the group of points clusters together close to the origin 
under their mutual gravitational attraction, producing a quite ir- 
regular mesh there. This effect is prevented when additional shape 
correction forces are invoked, as shown in the top panel. 



mass in the background can still be ignored in the evolution 
of the system. We calculate the gravity with the same tree 
algorithm used in the GADGET2 code, simply using the N- 
body approach with the masses and positions of the VPH 
particles. It turns out however that particles sometimes tend 
to pair up under their pairwise gravitational forces in the 
plain VPH scheme. This is related to the same defect dis- 
cussed in the context of Figure [2] If two particles are very 
close to a Voronoi wall, they can be moved still closer to- 
gether without increasing the hydrodynamic pressure force, 
so that a residual gravitational attraction (if not damped out 
by the gravitational softening) can move the particles very 
close together, with problematic consequences for the stabil- 
ity and accuracy of the scheme. A time-dependent gravita- 
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tional softening (|Price fc Mon aghan 20Q3), where the soft- 
ening is somehow tied to the size of the cell associated with 
a particle, may ease the problem, but is unlikely to cure 
it completely. However, the shape correction forces we in- 
troduced in Section [3.21 can nicely solve this problem. This 
effect is illustrated in Figure 1161 where we show the central 
mesh geometry for a two-dimensional version of the Evrard 
collapse problem. In the following we consider therefore a 
calculation of the 3D Evrard collapse with our usual choice 
of /3 = 1.2 and /9i = 0.1. 

In Figure [T71 we show radial profiles of gas density, ra- 
dial velocity and entropy for the Evrard collapse at time 
t = 0.8, calculated with 24464 particles inside the initial 
sphere. We compare the VPH result (shown by blue circles) 
to results obtained with SPH (red diamonds), and com- 
pare these to results of a ID high-resolution PPM calcu- 
lation of the problem (soli d black line) provided to us by 
ISteinmetz fc Mueller! l| 19931 ). We find that the collapse is 
essentially equally well described with VPH as with 
SPH, with a slight advantage for VPH, which more 
accurately resolves the central density and captures 
the shock more sharply. However, these differences 
lie well in the range of changes one obtains for dif- 
ferent viscosity prescriptions, and therefore do not 
seem to be particularly significant. We also note that 
the extra viscous forces needed to maintain a regular par- 
ticle distribution in VPH do not introduce any unphysical 
features in the solution. In particular, the radial profile of 
the specific entropy shows no signs of extra cooling or heat- 
ing. 



6 CONCLUSIONS 

We have discussed a new fluid particle model where the den- 
sity is estimated with the Voronoi tessellation generated by 
the particle positions. Unlike in SPH, there is an auxiliary 
mesh, which adds complexity to the scheme. However, the 
use of this fully adaptive mesh offers a number of advantages. 
It offers higher resolution for a given number of particles, 
since fluid features are not inherently smoothed as in SPH. 
In fact, the tessellation techniques are probably close to an 
optimum exploitation of the density information co ntained 
in the particle distribution (|Pelupessv et al.l I2003T ). When 
even higher resolution is needed the method could be ex- 
tended by an adaptive particle refinement or splitting tec h- 
nique, for example similar as sug gested in (|Sprin gcl 2009). 

As a result, contact discontinuities can be resolved with 
one cell, and surface tension effects present in standard SPH 
across contacts with large density jumps are eliminated. This 
has further implications for the growth rate of fluid instabil- 
ities in inviscid gases. Furthermore, the free parameters in 
the density estimate of SPH, involving both the number of 
neighbours as well as the kernel shape, are eliminated, which 
can be viewed as a good thing since the optimum values for 
them are not known, and incorrect choices can invoke the 
well-known clumping instability in SPH. 

One somewhat problematic aspect of Voronoi-based 
particle hydrodynamics is that the noise in the scheme is 
quite sensitive to the level of mesh regularity. Flows with a 
lot of shear can readily develop Voronoi meshes with points 
that lie close to the surfaces of their Voronoi cells. In this 
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Figure 17. Evrard collapse at time t = 0.8 simulated with SPH 
(red diamonds) and the Voronoi scheme with shape correction 
forces (blue circles). The black line shows the results of a ID 
hig h-resolution PPM cal culation of the problem provided to us 
by ISteinmetz fe Muellerl jl993h . From top to bottom, we show 
radially averaged profiles of gas density, radial gas velocity, spe- 
cific entropy, and internal energy. 
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Figure 18. Energy evolution for the "Evrard collapse" simulated 
with SPH (black dashed line) and the Voronoi scheme with shape 
correction forces (red solid line). The small fluctuations in the 
total energy arise primarily as a result of the finite accuracy of 
the tree code used to calculate self-gravity, and are of similar 
magnitude for both cases. 



case the noise in the gradient estimates increases, and, more 
importantly, it becomes difficult to safely prevent particle 
interpenetration, since closeness to a wall of the tessellation 
always implies that there is a second point on the other side 
of the wall which is also close, i.e. in other words, that a 
close particle pair is present. 

Higher order density estimates might solve this prob- 
lem, but they would have to be introduced already into the 
Lagrangian, leading to much more complicated equations of 
motion that may be intractable. We therefore explored two 
different approaches for keeping the mesh relatively regular. 
One is simply based on trying to formulate additional ar- 
tificial viscosity terms such that the viscosity tries to make 
cells 'rounder'. Whereas this shows some success, it does not 
succeed in all situations, particularly in strong shear flows 
where the artificial viscosity needs to be very low. A more 
radical approach also explored is to add correction terms 
to the underlying fluid Lagrangian with the aim to penalize 
strong deviations from regular mesh geometries. Our goal 
was to impose small, non-dissipative correction forces that 
maintain a proper mesh geometry. Thanks to the Lagrangian 
formulation, the required form of the correction forces to re- 
tain fully conservative behaviour can readily be derived, and 
the fluid motion under these forces shows the desired proper- 
ties. However, if the correction terms are too large, one risks 
deviations from the proper hydrodynamic solution. Further 



experimentation will be required to identify the optimum 
setting of these parameters. 

The Voronoi-based fluid particle approach can be rel- 
atively seamlessly integrated into an existing SPH code, 
provided a tessellation engine can be added in an appro- 
priate fashion. Other aspects of the physics (in particu- 
lar self-gravity, an additional collisionless component, radia- 
tive cooling, star formation, and feedback processes) can be 
treated in essentially identical ways as in SPH. This makes 
it possible to readily apply Voronoi particle hydrodynam- 
ics to problems of interest in cosmological structure forma- 
tion. In general, our first results suggest that VPH is su- 
perior to SPH, albeit at much increased complexity. How- 
ever, it is at present still unclear whether it is competitive 
with finite volume hydrodynamics carried out on a similarly 
constructed Vo ronoi mesh, as realized in the AREPO code 
l|Springelll2009l ). To elucidate this point further, we are in the 
process of carrying out galaxy collision simulations as well as 
cosmological structure formation simulations with our new 
technique and will report the results in forthcoming work. 
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Figure Al. Sketch of the integration variables 
APPENDIX A: DIFFERENTIAL OPERATORS ON VORONOI MESHES 

In this Appendix, we collect some useful formulae for discretized versions of differential operators on Voronoi meshes, such as 
the gradient or the divergence, and we test the accuracy of the gradient estimate as a function of the regularity of the Voronoi 
mesh. 



Al Gradient 



The cell averaged gradient of any quantity can be estimated via Gauss' theorem. One can use Gauss' theorem on a product 
of a constant vector times a scalar field and arrives at: 



V 



— / Vd>dV = — 



v 



V 



>dS, 



(Al) 



which can be used as one way to derive an estimate of the local gradient by approximating the value of cf> on the surface of a 
cell with the arithmetic mean between the cell and its neighbours. However, one can also circumvent the problem of finding 
a proper value for <j> on the surface by using a different starting point. We now use Gauss' theorem on (1 ■ r) V</>, where 1 is 
the unit vector. This yields: 

(A2) 



Vcj>dV = I r{V<j)-&S)- J rAcj>dV 

dV 



With the help of Rij — \rj — n\ as the distance between points i and j, and r x as r projected onto the plane of the face Aij 
(for definition see Fig. Q]and Fig |Al|) . r becomes r — r x + (n — Vj)/2 + TVj, and then the right hand side can be discretized 
for the Voronoi mesh as follows: 

(A3) 

rA<j)dV (A4) 

rA^dV (A5) 

(A6) 
(A7) 

The second term vanishes for linear scalar fields. It is therefore only a second order correction that becomes negligible for 
sufficiently smooth fields if the points lie near the centroids of their cells so that J (r — r;) dV . 

j>j — cj>i), so that we obtain for the gradient estimate (V I ) i at point i 



/ r(V<t>-dS) = / r (' S7( t ) - ds ) - / rA<l>dV 

■lav jjti J A i:j Jvt 
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We note that application of this gradient estimate to the Euler momentum equation in the form 

nrnri — —Vi VPi 

yields 



•2> 



{Pi-Pi)??- + (Pi+Pi, 



(A10) 



(All) 



which is consistent with the expression derived directly from the Lagrangian. 



A2 Divergence and curl 

To estimate the divergence and curl of the velocity which we need for the viscosity calculation of Section [2]4] we use the same 
reasoning as in[ 



(rV (V«) + Vv) dV = / r(Vw-dS). 

Vi J dV 



Provided that ^ J v V x (V x v) vanishes for a linear field we define our estimator for the divergence operator as 



ji 1 * 



I 1 Cii 



Similarly, the curl estimator can be defined in the form 
(VxB) i = -i^ 



iB,-B,) (^ lJ + gL 



(A12) 



(A13) 



(A14) 



where B denotes some vector field. Again applying Gauss' theorem to V(j> the Laplacian of a scalar function <f> can be computed 
as 



A(j>dV = ^2 (V<£-dS) A 
3 yi ■' A 



Ri 



(A15) 



A3 Accuracy of the gradient 

To test the accuracy of the numerical gradient estimate, we assume a quadratic model function 
and Hesse matrix, of the form 

1 T 

4>(r) = 4>o + At + —r Br. 

For definiteness, we set B — 67, where I is the identity matrix. We then populate a box of unit length on a side with a set 
of points, and evaluate the function 4>(r) at the coordinates of each of the points. After constructing the Voronoi tessellation 
for the point set, we then estimate the local gradient for each cell based on 



(r) with constant gradient 



(A16) 



(V0)< = y_ ^{4>j - <t>i) Aij ( ie tj + 



(A17) 



and alternatively also based on a simpler version of this formula where the terms proportional to dj are omitted, which 
corresponds to the simplest version of a Green-Gauss gradient estimate. We use three different point distributions with 4096 
points in a box of size unity. First we use a (i) random Poisson point distribution, a (ii) relaxed point distribution obtained 
from the VPH scheme where each cell has the same volume (obtained from the top of Fig. [6|, and (iii) a distribution relaxed 
with PPO (obtained from the bottom of Fig. [6} where in addition very round cells were produced in which the majority of the 
points lies close to the geometric centres of the cells. In all three cases we compare the magnitude of the estimated gradient 
vector to the magnitude of A, and we plot the median of the relative error as a function of 6/|A|. To exclude boundary effects, 
only cells whose neighbours do not overlap with the box boundary are considered in the measurement. 

In Figure IA21 we show the results. The panel on the left gives our adopted gradient estimate, while the panel on the 
right is for the simpler version of the Green-Gauss gradient estimate. Interestingly, for 6 = 0, the error vanishes exactly, 
independent of the regularity of the Voronoi mesh. However, once the second order term starts to influence the measurement, 
i.e. for large values if b/\ A\, the more regular meshes clearly yield a lower error, as expected. In all cases, the gradient estimate 
that includes the Cij term is superior to the simple Green-Gauss gradient estimate. In particular, only when it is included, a 
vanishing error for a linearly varying field is obtained. 
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Figure A2. Median relative error in the gradient estimate obtained either with our default formula (left) or with the simpler Green- 
Gauss estimate where the terms are omitted (right). We show results for three different types of point sets, a Poisson sample (black), 
a regularized distribution where each Voronoi cell has equal volume (red), and a regularized distribution where in addition the cells are 
quite 'round' and regular (blue). The accuracy is measured as function of the strength of a second order variation in the underlying field. 



APPENDIX B: CONTROLLING THE SHAPE OF CELLS 



As outlined in Section I3.2I we modify the fluid Lagrangian slightly to include factors that penalize highly distorted cell 
shapes. If such shapes occur, we want small adjustment forces to appear that tend to make the mesh more regular again. 
These adjustment forces need to preserve energy and momentum conservation of the scheme, which will automatically be 
the case if they are derived from a suitably defined Lagrangian or Hamiltonian. In this Appendix, we derive the equations of 
motion for the Lagrangian 
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m k r\ 



P k V k 
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(Bl) 



where the factor in square brackets disfavours displacements of points from the geometric centres of their cells, and the factor 
in curly brackets disfavours cells with large aspect ratios. 
We define the centroid of a cell as 



_L_ 

Vk 



rxk(r)dr, 



where \k is the characteristic function of cell k. The shape of a cell is measured via the second moment 
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((r - s k ) 2 ) k = — J (r- s k ) 2 X k(r) dr 



(B2) 



(B3) 



Here d counts the number of dimensions, i.e. d = 2 for 2D and d = 3 for 3D. The factor v£ is hence proportional to the 
'radius' R k = Vu of a cell squared, fie measures the strength of the effect of displacements of points from the centroid of 
a cell, while /3i is the corresponding factor for the aspect-ratio factor. The constant /?2 is only introduced to prevent that 
even round cells lead to a significant enhancement of the thermal energy. For perfectly round cells, we expect in 2D roughly 
circles for which w k = V 2 ' d /(2Tv), hence we pick P2 = l/(27r). In 3D, we have approximately spheres instead and we pick 
/3 2 =3/5(3/4tt) 2/3 . 

Note that this Lagrangian is only a function of the point coordinates for given entropies, so the equations of motion for 
conservative dynamics are perfectly well defined, even though they lead to more lengthy expressions than in the standard 
case. We first obtain the following Lagrangian equation of motion: 
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(B4) 



where the empty square and curly brackets are notational short-cuts for the corresponding terms in the original Lagrangian. 
We note that we can use the identity 



j-PkVk = 
dri 



vp dV k 



(B5) 
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We also note that in the second and third terms of equation (|B4[) we encounter partial derivatives of Vk, which we can combine 
with the first term into a more compact form. This allows us to write the equation of motion in the form 

dri 



dr 

k k 

where we have defined 
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We already know an explicit expression for dVk/dri, but we still need to derive such a thing for the derivatives of (rk — s k ) 2 
and w k . Let us first deal with the term involving Qk in the equations of motion, i.e. 



(rmri) Q = 



k ' k 



d{r k - s k ) 
dri 



(r k - s k ). 



(BIO) 



Here the exponent T stands for the transpose, and the notation ^ is the Jacobian matrix with elements = -§p~- 

Based on the definition of s k in terms of the characteristic function we find 
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For the derivative of the characteristic function we can use a result from ISerrano fc Espanoll l|200lf ) and write 
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which is based on approximating the characteristic function with 
Xk(r) = 
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which becomes exact in the limit a — > 0. Putting these results into equation (|B10|) one gets 
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We can now identify the area of a face between two cells as 



Aij — Rij 



dr 



XiXi 



and the centroid of the face as 
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i XiXi 
dr — Tr-r. 



Furthermore, we define a second-order tensor of the face relative its centroid as 
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With these definitions, we can rewrite equation (|B14|) as 
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where we introduced the further short-cut 
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(B19) 



We note that the second term in this equation can be absorbed in yet a further redefinition of P£, which we will exploit later 
on. We next consider the term in the full equation of motion that involves the L k factor. This is given by 
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We now define a further moment for each cell face, namely the vector-valued quantity 
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Note that g ij always vanishes in 2D but can be non-zero in 3D. With this definition, we can rewrite equation (|B20[) as 
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where we have defined the short-cut 
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Again, the first term involving dV k /dri can be absorbed into a redefinition of P k . Putting everything together, the complete 
equation of motion can then be written as 
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V k V k 

While a bit lengthy, this can be straightforwardly calculated for the VPH scheme. Nevertheless, we want to add a brief note 
on how to compute the Tensors Ty , which is done as part of the mesh construction. We have 

<(r - sf) k = ((r - r ) 2 ) fc - (r„ - s) 2 (B26) 
for any reference point ro- Suppose we have a triangle in 2D given by (ro, ri, r 2 ), then the moment can be obtained as 

((r - r o) 2 ) k = g [(ri - r )(ri - r ) T + (ri - r )(r 2 - r ) T + (r 2 - r )(r 2 - t- ) t ] . (B27) 
Similar relations hold for 3D and can be exploited for an efficient calculation of the tensors Ty and the vectors . 



